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1 INTRODUCTION 



Over the past decade, the cold dark matter model, comple- 
merrted with dark energy (ACDM), has established itself as 
the theoretical framework of choice to describe the formation 
of structures in the Universe. Whilst dark matter dominates 
the dynamics of structure formation on large scales, a host 
of observational evidence indicates that the situation is re- 
versed on galactic scales dKassin. de Jong. Weiner 2006l l . In 
particular, the inner structure of galactic d ark matter ha - 
los must be affected by baryonic processes ijMerritt 20061 ). 
Whereas it is certainly true that critical aspects of the 
baryonic physics of galaxy formation are still poorly un- 
derstood and beyond the reac h of direct numerical sim- 
ulations ijSpringel et al. 20051 1. the overwhelming major- 
ity of observational constraints unfortunately comes from 
the electromagnetic emission of these baryons. Therefore, 
significant progress in our understanding of the theory of 
galaxy formation and evolution depends on our ability to 
perform hydrodynamical cosmological simulations. At first, 
these will necessarily be complemented by appropriate sub- 
grid modelling, presumably based on knowledge gained from 
hydrodynamical simulations performed on smaller scales 



ABSTRACT 

The largest galaxies acquire their mass early on, when the Universe is still youthful. 
Cold streams violently feed these young galaxies a vast amount of fresh gas, resulting 
in very efHcient star formation. Using a well resolved hydrodynamical simulation of 
galaxy formation, we demonstrate that these mammoth galaxies are already in place 
a couple of billion years after the Big Bang. Contrary to local starforming galaxies, 
where dust re-emits a large part of the stellar ultraviolet (UV) light at infrared and 
sub-millimetre wavelengths, our self-consistent modelling of dust extinction predicts 
that a substantial fraction of UV photons should escape from primordial galaxies. 
Such a model allows us to compute reliably the number of high redshift objects as a 
function of luminosity, and yields galaxies whose UV luminosities closely match those 
measured in the deepest observational surveys available. This agreement is remark- 
ably good considering our admittedly still simple modelling of the interstellar medium 
(ISM) physics. The luminosity functions (LF) of virtual UV luminous galaxies coincide 
with the existing data over the whole redshift range from 4 to 7, provided cosmolog- 
ical parameters are set to their currently favoured values. Despite their considerable 
emission at short wavelengths, we anticipate that the counterparts of the brightest 
UV galaxies will be detected by future sub-millimetre facilities like ALMA. 

ijSlvz et al. 20051 1. Indeed, we are seeking to describe highly 
nonlinear mechanisms that are coupled together, such as ra- 
diative gas cooling, star formation and feedback. Such com- 
plexity makes it, for instance, very unlikely that even when 
armed with a satisfying model for supernovae explosions, 
one can describe analytically their interplay with a clumpy, 
multiphase ISM and/or intergalactic medium (IGM). 



Perhaps the most significant fiaw of all galaxy formation 
models which rely on the results of pure dark matter simula- 
tions is their neglect of the filamentary cold gas flows which 
feed galaxies ijPekel. Birnboim 2006l : lDekel et al. 2009l l. Gas 
accreted in this mode substantially dominates the accretion 
budget. Indeed, cold flows are the unique mode of gas accre- 
tion for galaxies hosted by all dark matter halos with mass 
less than 5 10^^ Mq (i.e. the vast majority of galaxy size 
halos). They are reported in every high resolution cosmolog- 
ical hydrodynam ics simulation of galaxy formation w hich 
includes cooling ijKeres et al. 20051 : lOcvirk et al. 2o"08l l. re- 
gardless of the technique employed to solve the Euler equa- 
tions. Such an oversight is predicted to have drastic conse- 
quences on the ability of these models to match star forma- 
tion rates and supernova driven gas outflows. Arguably, the 
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Figure 1. A global composite multiscale view of the Mare Nos- 
trum simulation described in the text. The blue colour traces the 
gas temperature, green represents the gas density and red shows 
the dark matter density. The box size is 50h~^ co-Mpc on a side 
for the top panel and 10h~^ comoving Mpc for the bottom one. 
Galaxies appear as pinkish yellow smallscale disk-like structures. 



most straightforward observational tracer of the cold flow 
scenario is the UV flux emitted by galaxies, as it mainly 
originates from newly born stars. However, the matter is 
rapidly complicated by the presence of dust, which is gen- 
erated on supernovae timescales (« 50 million years) and 
is very efficient at absorbing UV radiation. The spatial dis- 
tribution of dust grains relative to stars is therefore a cru- 
cial ingredient in any estimate of the UV luminosity func- 
tion a nd it is notoriously difficult to predict (semi -) analyti- 
cally ijPevriendt et al. 19991 : iFontanot et a?. 2009l l. Accord- 
ingly, we address this key issue numerically. 



2 SETUP AND EXTINCTION MODEL 

The evolution of a cubic cosmological volume of 50h~^ Mpc 
on a side (comoving), containing 1024^ dark matter par- 
ticles and an Eulerian root grid of 1024^ gas cells is fol- 
lowed. A ACDM concordance universe {yi,n ~ 0.3, = 
0.7, /i = Ho/[lQQ kms-^Mpc-i] = 0. 7, as = 0-9, n = 1, i.e . 
the WMAP 1 year best fit cosmology ijSpergel et al. 2003h . 
is adopted, resulting in a dark matter particle mass mp = 
1.41 X IO^Mq an d we use an a daptive mesh reflnement 
(AMR) technique ijTevssier 2002l l to keep our spatial res- 
olution flxed at around 1 h~^kpc in physical coordinates. 

The matter density flelds are evolved from z = 120 
to z = 4, and we output about 34 snapshots regularly 
spaced in the expansion factor. In each snapshot, we iden- 
tify dark matter halos as well as the subhalos they con - 
tain using the AdaptaHOP algorithm dAubert et al. 20o3 ). 
and only keep groups with more than 100 particles. We 
measure the basic physical properties of these groups and 
subgroups, such as their masses Mhop and radii iihop. Key 
physical processes for galaxy formation such as radiative 
cooling, UV background radiation, star formation and su- 
pernovae fe edback are also implemented self-consistentl y 
in RAMSES ijRasera. Teyssier 20061 : iDubois. Teyssier 20081 ). 
Gas cools down to a minimum of 10* K at a rate which de- 
pends on the local metallicity in each grid cell, and turns into 
stars on a constant timescale off* = 2Gyr. Stars then evolve, 
releasing an amount of metals and energy into the ISM spec- 
ifled by a Salpete r initial mass function and the physics of a 
Sedov explosion ijDubois. Teyssier 2008l l. In some cases, the 
combined supernovae blow part of the ISM back into the 
halo's hot phase and further into the IGM. For each time 
snapshot of the simulation we associate the stars located 
within each dark matter substructure detected by Adap- 
taHOP to a different galaxy. Note that, proceeding this 
way, we account for all the stars created in the simulation. 
By keeping track of the stellar content of each galaxy, as 
a function of age and metallicity, and knowing the galax- 
ies' gas content and chemical composition, we then compute 
the (possibly extinct) spectral energy distribution of each 
galaxy, along the line s described in the Stardust model 
(Devriendt et al. 1999l l. 

More speciflcally, we calculate the spectrophotomet- 
ric properties of galaxies including internal extinction of 
starlight by dust. We model dust absorption within each 
galaxy using the empirical calibration of Guiderdoni & 
Rocca-Volmerange (1987), in which the optical depth scales 
with metallicity and in proportion to the column density of 
hydrogen along the line of sight. The optical depth, r, is 
therefore given by: 

where A is the wavelength, and {Ax/ A\)zq is the extinc- 
tion curve for sol ar metallicity, which we take to be that 
of the Milky Way ijCardelli et al. 1989l l. We adopt a scaling 
with m etallicity of s=1.6 for A ^ 2000A and s=1.35 for A 
2000 A ijGuiderdoni. Rocca-Volmerange 1987l l. We calculate 
the extinction for subsamples of galaxies in logarithmic bins 
of stellar mass (100 objects per bin), based on the actual 
distribution of gas and metals in each galaxy and averaging 
over lines of sight. The resulting distributions are then used 
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to draw random extinctions for galaxies with similar stellar 
masses. The total number of galaxies varies between 43000 
and 86000 between z=7 and z=4. 



3 THE UV LF OF HIGH-Z GALAXIES 

Observationally, the largest samples of UV emitting galaxies 
at « > 3 are selected according to a colour-colour criterion 
that maps the Lyman break feature in their spectral energy 
distri butions (hence their name "Lyman-Break Galaxies" or 
LBG ijSteidel ei ai. 1999l ll. This now well-established tech- 
nique allows the detection of large numbers of galaxies per 
field of view using relatively little obs erving time, and does 
not strongly bias the selection process ijBouwens et al. 20061 : 
iBunker et al 20061 : 1 Ilbert et al. 2 005^. Such data sets pro- 
vide us with the best sample of high-z galaxies against which 
to perform a meaningful statistical comparison. The most 
robust and simple statistics to consider for such a compar- 
ison is the LF of galaxies. This is the main focus of this 
letter. 

Like the observations, which are limited by the resolving 
power of the telescope/instrument combination with which 
they are obtained, the LFs measured in numerical simula- 
tions suffer from finite box size and resolution effects, so 
their domain of validity is restricted to a certain range of 
magnitudes. In the simulation presented here (Fig 1), dark 
matter halos smaller than ^ IO'^Mq in mass are barely re- 
solved in terms of the hydrodynamics of their gas content, 
so galaxies hosted by smaller halos are absent, which in- 
troduces a non-physical cutoff at the faint end of the LF 
(Fig. 3). As spectral energy distributions of galaxies are com- 
puted as the weighted sum of single stellar populations, we 
avoid being dominated by Poisson noise by restricting our 
results to galaxies containing at least 10 star particles. At the 
bright end, objects become rarer and rarer, until we reach 
the point where we run out of bright galaxies. This effect is 
somewhat amplified by our finite simulated volume, which 
is small compared to the size of the current observable Uni- 
verse: at some point, the statistical error on the LF becomes 
of the order of the measurement itself. Vertical and horizon- 
tal dashed lines on Fig. 3 mark these numerical limitations 
respectively, clearly delineating the region where our LFs can 
be trusted. Note that, although approximate, these are con- 
servative boundaries. A truly accurate determination would 
require running several simulations encompassing larger vol- 
umes and featuring a higher mass resolution. Finally, despite 
our resolution, still short of being able to provide a fully 
self-consistent model for re-ionization by about an order o f 
magnitude in mass ijGreff et al. 20081 : IWhalen et al. 20041 ). 
the UV backgroun d radiation ha s little effect on the galax- 
ies discussed here ijGnedin 2000l l since the minimal circular 
velocity of their host d ark matter halo is larger than 25 km/s 
ijOkamoto et al. 2008l l. 

Setting aside the star formation rate and stellar initial 
mass function, the key issue to derive correct UV LFs lies in 
our ability to estimate dust extinction accurately, since dust 
grains very easily absorb light in this wavelength range. Be- 
cause a full modelling of the formation and evolution of the 
dust grain population is beyond the reach of current simu- 
lations, we have tried a couple of phenomenological models 
to assess the impact of extinction on our results. First, as 



a sanity check, we use the most widespread method in the 
literature, which assumes a uni versal shape for the extinc- 
tion curve ijCalzetti et al. 1993 l. and then fits its normal- 
ization, E(B— V), to obtain the bes t possible match to the 
observational luminosity function ijNight et al 2006l l. This 
means that a unique "average" extinction is assigned to ev- 
ery galaxy for which we compute a UV magnitude. We be- 
lieve this is, at best, a very awkward state of affairs, since it 
is empirically known that, at least in the local Universe, the 
value of the extinction is correlated to the UV magnitude, 
such that brighter ga laxies have larger values of E(B — V) 
ijHeckman et al. 2005l l. However, as pointed out by Night et 
al. (2006), invoking luminosity dependent extinction can be 
seen as a cheat, as it is equivalent to hiding the discrep- 
ancy between simulated and observed LFs in the extinction 
law. As a more ambitious alternative, we therefore measure 
the distribution of gas and metals on the line of sight to- 
wards each of the star particles in our simulated galaxies, 
which allows us to calculate the optical depth as a function 
of wavelength for each of these lines of sight (see Fig. 2). 
The main caveat of this second approach, is the assump- 
tion that dust properties only depend on these two quanti- 
ties (at least for absorption of UV lig ht), but it was shown 
ijGuiderdoni. Rocca-Volmerange 19871 1 that this seems to be 
a fair description for extinction in local galaxies. Arguably 
the most important feature of such a dust model is that it 
provides a self-consistent method to compute the extinction 
law directly from the simulation itself. 

The two dust models produce LFs that match the 
observed one quite well within the whole range of mag- 
nitudes it covers, as observational error bars plotted 
on fig. 3 only roug hly estimate cosmic variance errors 
ijBouwens et al. 2008l l. We discuss our results in terms of 
the most commonly used analytical parameterisation of 
the galaxy LF, i.e. the Schechter function, examples of 
which are plotted on fig. 3. This description relies on three 
parameters, M* and a, i.e. the normalisation, char- 
acteristic magnitude and faint-end slope of the LF, re- 
spectively. The preferred Schechter fit for our LF at z=4 
(i.e. the green curve on the bottom right panel of fig. 3) 
yields a faint end slope of 1.6 and a characteristic mag- 
nitude of 20.9, comp ared to 1.73+/0.05, 20.98+/0.1 for 
the data respectively ijBouwens et al. 2007f ). Our normali- 
sation of 0.003 Mpc"^ appears a bit high with respect to 
that measured in the data (0.0013+/0.0002 Mpc^). How- 
ever, this is the parameter of the Schechter function which 
is the most sensitive to a change in cosmological model. 
In other words, the error we make when we rescale our 
LFs from the Wi lkinson Microwave Anisotropics Probe 
(WMAP) 1 year ijSpergel et al.. 2003l l to the WMAP 5 
year parameters (green curves on fig. 3, flm ~ 0.26, Qa = 
0.74, ft = Ho/\100kms ~^Mpc~^] = 0.72, crs = 0.8, n = 0.96) 
ijPunklev et al.. 2009l l mostly affects so that the factor 
of about two difference is not statistically discrepant. To be 
more precise, in order to recast the results of our simula- 
tion for a different set of cosmological parameters, we have 
used the Press & Schechter (1974) formalism to calculate 
the dark matter halo mass function at a given z in both cos- 
mologies. Assuming an occupation number of one galaxy per 
dark matter halo, we then divide our luminosity functions 
by the ratio of the dark matter halo mass functions. Given 
the mass resolution of our simulation, we expect this correc- 
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Figure 2. A typical field of view of 1 comoving Mpc on a side in 
the Mare Nostrum simulation shown in composite rest frame true 
colours (R, B, and V filters). Colours of the stellar population 
are shown in the left panel without dust extinction, and with 
extinction in the right panel. Note the presence of dust lanes in 
the two largest galaxies, testimony of the highly non- homogeneous 
distribution of dust and stars in the ISM. 

tion to be quite accurate at higher redshift (above z=6) and 
to slowly degrade as more massive haloes begin to form and 
start hosting a larger number of galaxies on average. This 
explains, at least partially, the larger discrepancy (at the 
~ 2 sigma level once cosmic variance is taken into account) 
between model and observed LFs at intermediate and low lu- 
minosities for z < 6. Note that these considerations do not 
apply to the brightest objects, which are much more scarce 
and have a much higher probability of being unique in their 
host dark matter halo. As a result, the correction remains 



valid for these objects even for z<6, as can be seen on fig. 
3. 

In a nutshell, we find that our simulated LFs at z~4 are 
in respectable agreement with the available UV data and 
robustly predict a faint-end slope with a value a~ — 1.6 
compatible with the findings of Steid el and collaborators 
ijSteidel et al. 19991 : iReddv et al. 200I ) ( q » -1.6) but not 
with the slope found in the Subaru data ijOuchi et al. 20[i3 ) 
(q ^ —2.2). Our LFs also show quite an important amount of 
evolution, as measured by M*, which varies from 19.8 at z=7 
to 20.9 at z=4. This strong evolution is a direct consequence 
of the rapid decrease of the cosmic star formation rate den- 
sity from z = 4 to 7 in our simulat ion and is, to a large ex tent, 
unaffected by resolution effects ijRasera. Teyssier 2006h . At 
z>6, recent analyses of the Hubble Space Telescope Ultra 
Deep Fie ld and the Great Ob servatories Origins Deep Sur- 
vey data ijBouwens et al. 20071 : B ouwens et a l. 2Qo3) suggest 
a slightly steeper faintend slope of a w —1.74, again in good 
agreement with our simulation as we find that the faint-end 
slope of our LF does steepen slowly with redshift. This is of 
key importance since ionising photons from low luminosity 
galaxies could possibly suffice to carry out the reionisation 
of the Universe, provided these latter are numerous enough. 
Reversing the argument and marching down in redshift, this 
raises the question of what mechanism is responsible for such 
a flattening of the faint end slope of the UV LF, which ex- 
tends down to the local Universe, since a shallow slope (a 
pa 1.22) is observed locally for z < 0.1 in the GALEX far- 
UV data ijWvder et al. 2005h . Supernova feedback could 
potentially drive powerful galactic outflows and quench star 
formation, but SPH simulations that explicitly include winds 
efficient enough to drive gas out of low mass halos, system- 
atically m^asuje_a_steeper faint end slope than we do (a ~ 
-2.0) ijNight et al 20061 ). Finally, at the bright end of the 
LF, our two dust extinction models yield different results. 
The self-consistent extinction model has a more pronounced 
cut-off, which better matches the shape of a Schechter func- 
tion. This feature stems from an intrinsically stronger dust 
attenuation that occurs in the most UV bright galaxies. We 
therefore predict that these very same UV beacons also are 
the most luminous IR sources in the sky. The forthcom- 
ing generation of sub-millimetre facilities, and in particular 
ALMA, will either validate or disprove this prediction. 



4 CONCLUSIONS AND DISCUSSION 

We have derived UV LFs of simulated high-z galaxies us- 
ing the largest hydrodynamical cosmological simulation per- 
formed so far. Including a self-consistent post-treatment of 
the extinction in our galaxies based on the abundance of gas 
and metals measured in the simulation, we have compared 
these synthetic LFs to the observed ones. We have studied 
how sensitive these LFs are to both cosmological parame- 
ters and dust absorption, using th e most commonly u sed 
extinction model in the community ijCalzetti et al. 19941 ) as 
a reference. Our self-consistent treatment of dust absorp- 
tion enables us to match the observed UV LFs. This match 
quantitatively compares to the fit obtained with the sim- 
ple model where all the galaxies are dimmed in the same 
way, provided that we set the cosmological parameters (in 
particular the normalization of the power spectrum and 
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Figure 3. Redshift evolution of tiie UV (1400 A) rest frame galaxy luminosity function measured in the Mare Nostrum simulation. 
Vertical dashed lines indicate mass resolution, horizontal dashed lines volume resolution. Purple diamonds indicate nonextinguished 
WMAP-1 cosmology LPs, blue stars LPs extinguished with a dust content scaled using the measured a bundance of metals and gas in the 
simulation. Red diamond LPs stand for galaxies uniformly extinguished with a l|Calzetti et al. 1994 ) law in a WMAP-1 universe, and 
green stars for LPs extinguished as the blue star ones but res eated to the currently favoured WMAP-5 cosmology. Black triangles and 
solid lines mark data gathered by Bouwens and collaborators ijBouwens et al. 20081 : iBouwens et al. 2007^ 



its rollinR index) to their best fit WMAP 5 years values 
ijPunklev et al.. 2009l l. However, contrarily to the Calzetti 
extinction model, such an agreement is achieved without the 
introduction of an extra free parameter, i.e. the average ex- 
tinction of the galaxy population as a function of redshift can 
be retrieved from the simulation itself. Moreover, this self- 
consistent model receives support from (i) the extrapolation 
of low z observa tions which show that extinction depends on 
UV luminosity ijHeckman et al. 2005l l (ii) the mounting evi- 
dence for an extreme ly low metallicity in faint z~7 galaxies 
ijBouwens et al. 2009l l. which both favour a non- uniform de- 
gree of extinction throughout the magnitude range spanned 
by the LF. 

Although current observational surveys still probe small 
volumes of the high redshift universe, and are therefore 
prone to suffer from cosmic variance at the bright end of 
the LF, they seem to go deep enough to yield reliable esti- 
mates of its faint end slope in the UV, up to z~7. At fixed 
cosmological parameters, the key to matching UV luminosi- 
ties with any model of galaxy formation is the amount of 



dust extinction present in galaxies. However, the absorbed 
UV light is reprocessed by dust, and astronomers should 
be able to detect the far infrared counterparts of low lu- 
minosity primordial galaxies, provided extinction strongly 
affects galaxies across their entire luminosity range. Future 
wider surveys together with the advent of infrared and deep 
millimetre observations from the JWST and ALMA respec- 
tively should therefore give us a robust handle on where 
energy from star formation is really coming out at these red- 
shifts. These measurements will, in turn, confirm or inval- 
idate our improved treatment of dust attenuation. Indeed, 
we predict that this attenuation should strongly correlate 
with UV luminosity, with low luminosity galaxies almost 
dust free at z~7, and massive galaxies extinguished by 1.5 
to 2 magnitudes. Hence, only the tip of the iceberg of pri- 
mordial galaxy formation should be visible at rest-frame in- 
frared wavelengths if our model is correct: the bulk of the 
galaxy population will mostly be UV bright with a few iso- 
lated monsters shining across the entire wavelength range. 
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